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1.   INTRODUCTION 

Considerable  literature  exists  describing  the  static  characteristics 
of  aircraft  near,  or  at,  stall,  including  detailed  analyses  of  stall 
tendencies.   However,  there  appears  to  be  little  information  on,  or 
analytic  treatment  of,  aircraft  dynamics  in  the  region  around  stall,  and 
especially  post-stall  time  histories  in  which  divergent  oscillations 
arise.   The  latter  phenomenon  is  quite  interesting,  since  experienced 
pilots  apparently  are  familiar  with  the  "rocking  chair"  or  "porpoising" 
stall,  where  the  aircraft,  while  under  constant  elevator  control,  alter- 
nately enters  stall,  then  recovers,  leading  to  a  vertical  oscillation  of 
the  aircraft  nose  about  some  point. 

Normally  we  define  the  airplane  stall  angle  as  that  angle  where  the 
vehicle's  trimmed  lift  curve  (C  vs  a)  has  a  global  maximum.   It  is  well 
known  that  the  lift  curve  then  exhibits  a  negative  slope  for  some  range 
of   a  above  the  stall  angle.   Therefore,  we  shall  also  use  here  a  ,  the 
stall  break  angle,  defined  as  that  angle  where  the  lift  curve  has  maximum 
negative  slope.   There  is  a  paucity  of  information  on  the  nature  of  the  lift 
curve  in  the  post-stall  or  stall  break  region,  since  experimental  wind 
tunnel  data  has  usually  been  smoothed  through  this  region  based  on  engineer- 
ing intuition.   In  addition,  an  actual  treatment  to  define  the  stall  break 
for  a  physical  aircraft  is  difficult  since  model  data  are  subject  to 
scaling  difficulties  and  aircraft  parameter  in  the  non-linear  region  is 
limited  at  present.  (Figure  1) 
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In  his  thesis,  Frederiksen  [1]  considered  the  problem  of  divergent 
oscillations  and  limit  cycles  in  this  region.   His  linear  analysis  was 
able  to  show,  numerically  for  a  single  aircraft,  that  a  sufficiently 
large  negative  slope  in  C    in  this  region  led  to  roots  for  the 

LiOt 

characteristic  equation  which  had  positive  real  parts  (i.e.  exponentially 
growing).   This  behavior  is,  of  course,  transitory,  and  disappears  by 
the  time  the  aircraft  returns  to  a  configuration  where  the  lift  slope 
is  positive.   However,  it  did  provide  an  insight  into  the  mechanism 
that  might  cause  these  oscillations.   Furthermore,  he  was  able  to  simulate 
numerically  time  histories  in  which  a  "rocking  chair"  stall  did  occur. 
However,  due  to  the  numerical  nature  of  his  work  and  his  concentration 
on  a  single  aircraft,  his  results  are  not  readily  extendible,  or  suitable 
for  general  analysis. 

This  investigation  attempts  to  analytically  model  the  post-stall 
behavior  in  a  general  class  of  aircraft.   Our  main  hypothesis  is  that,  if 
the  slope  of  the  lift  curve  is  sufficiently  negative  in  the  region  beyond 
the  stall  break,  then  the  aircraft  will  behave  as  if  it  sees  a  sharp, 
essentially  instantaneous  loss  of  lift.   If  this  is  the  case,  then  the 
lift  curve  can  be  more  conveniently  modeled  in  this  region  by  the  intro- 
duction of  a  hypothetical  "relay"  acting  at  the  stall  break  angle.   We 
assume  that,  except  for  this  nonlinearity ,   C  ,  C   and  0   can  be 
adequately  approximated  by  linear  functions.   Those  assumptions,  together 
with  conditions  on  the  magnitude  of  various  coefficients,  lead  to  a  system 
of  coupled,  constant  coefficient,  first  order  differential  equations,  with 


a  forcing  term  which  has  a  jump  discontinuity  across  a  line  in  the  phase 
plane.   Analysis  of  this  system  yields  a  necessary  and  sufficient  condition 
for  the  existence  of  a  limit  cycle,  and  transcendental  equations  predicting 
the  period  of  this  cycle.   These  conditions  are  then  applied  to  the  case 
considered  by  Frederiksen,  and  it  is  shown  our  model  yields  excellent 
agreement  with  his  results. 
2.   DERIVATION  OF  THE  MODEL 

The  aircraft,  assumed  to  be  in  longitudinal  flight,  will  be  described 
by  means  of  the  standard  coordinate  system  shown  in  figure  2.   In  this 
system,  the  equations  of  motion  become 

{  T  +  qS  [C   sin  a  -  C   cos  a  ] } 

u  =  -  g  sin  9  -  w  6   ,      (1) 

m 


qS  [C   cos  a  +  C  sin  a] 

w  = +  g  cos  9  +  u  9   ,  (2) 

m 


*r  cm  •  (3) 

y 


and, 


{  -  u  sin  a  +  w  cos  a}  (A) 

a  =  


(The  terms  in  (1)  -  (4)  are  defined  in  the  Appendix.) 

These  equations  are  obviously  non-linear,  even  though  the  coefficients 
may  be  linear  or  constant.   In  addition,  though,  in  general,  over  wide 
ranges  of  variation,  the  coefficients  are  non-linear  functions  of  the 
variables.   In  most  applications  where  the  range  of  the  variables  is 
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Figure  3 


small,  they  may,  however,  usually  be  adequately  approximated  by  linear 
functions.   As  is  obvious  from  figure  1,  this  is  not  the  case  though 
for  C   in  the  neighborhood  of  the  stall  break. 

We  shall  proceed  under  the  assumption  that  any  motion  of  the  air- 
craft takes  place  in  a  small  region  about  the  stall  break  angle.   Further- 
more we  shall  assume  that  CM  and  C^  can  still  be  approximated  linearly 
in  this  region.   We  suppose  the  aircraft  is  trimmed  for  constant  attitude 
and  angle  of  attack,  denoted  by   6   and   a   respectively,  and  that  the 
elevator  angle  and  thrust,  denoted  by  6   and  T,   are  constant.   The 
velocity  of  the  aircraft  in  the  trim  configuration  is  denoted  by  v 
Frederiksen's  results  do  not  show  a  significant  variation  in  v  ,   compared 
to  that  in  a  and   6   during  limit  cycle  behavior,  so  we  assume 

v  =  v  . 
o 

Considering  now  perturbations  off  the  trim  configuration,  we  let 

Aa  =  a  -  a   ,   and  A 6  =  9  -  0 
o  o 

As  noted  above,  we  will  approximate  C   and  CL   linearly,  so  that 

CD  "  CD0  +  Sa  A°   >  (5> 


and, 


CM  =  "  CM0  "  CMa  Aa  "  2v~  CMe  A< 


~-  C...  Aa  -  CMA   6    .  (6) 

2v    Ma       M6   e 
o  e 


The  key  ingrediant  of  this  investigation  is  modeling  C   at  the  stall 


break.   Reference  to  figure  1  indicates  that,  if  the  negative  slope  in 
the  region  past  the  stall  break  is  sufficiently  steep,  then  an  aircraft 
entering  this  region  experiences  an  almost  instantaneous  loss  in  lift. 
This  is  essentially  equivalent  to  the  action  of  a  hypothetical  "relay" 
in  the  lift  curve  function  C   .   Thus,  we  shall  model  the  lift  curve  as: 

Li 


C 

JL  '    "1.0   "L6   ue   2v   ~L9 

o 


C,  =  C^  +  CT  «.   6   +  -7T-     CT  ;  Af 


+  CLa  Aa  "  CLA  Sgn(a  "  ac)   •  (7) 

c 

where  sgn(x)  =  |x|/x  ,  x  ^  0  ,   sgn(0)  =  0,   and  a   denotes  the  stall 
break  angle.   Figure  3  shows  a  sketch  of  a  C   curve  generated  by  (7) 

Li 

superimposed  on  the  C   curve  from  figure  1.   (Note  that,  to  facilitate 

Li 

analysis,  we  have  assigned  algebraic  signs  so  that,  in  normal  cases,  all 
the  coefficients  are  positive.) 

Note  now  that  u  and  w  which  appear  in  (4)  can  be  replaced  by 
their  equivalent  expression  in  terms  of   (1)   and   (2).   If  we  do  this, 
we  find  (since  u  =  v  cos  a,  and  w  =  v  sin  a): 

a  =  -  -  sin  a  -  ^  C,  +  K   cos(6  -  a)  +  6  (8) 

mv         mv  L   v 

But  now,  using  our  assumption  on  v,  and  the  fact  that 

•     .       .... 

a  =  Aa  ,   and  9  =  A9  ,  and  9  =  A9  , 

we  have, 

Ad  =  -  -^—   sin(a  +  Aa)  -  ^—   CT  +  &-   cos(9   -  a  +  A9  -  Aa) 
mv       o         mv   L   v       o    o 
o  o       o 

+  A9  ,  (9) 


and, 

Ae  -  ass.  c     .  (10) 

y 

But  then,  according  to  the  assumed  forms  we  have  for  C   and  C   , 

Li  ri 

we  see  that  Aa  and  A6   depend  only  on  Aa,  A0,  and   A6   (since 

(9)  can  be  used  to  eliminate  the  Aa   term  from  Cw  .)   Thus,  for  the 

M 

remainder  of  our  investigation,  we  shall  be  concerned  only  with 
equations  (6),  (7),  (9),  and  (10). 

Now,  under  the  assumption  that  Aa  and   A6   are  small,  we  see: 

sin(a  +  Aa)  =  sin  a  +  (cos  a  )  Aa   , 
o  o         o 


and, 


cos(6      -  a     +  A6   -  Aa)   =   cos(9      -  a   )    -  sin(6     -  a   s    AA 
o  o  o  o  o  o)    A8 


+sin(6      -a)Aa        , 
o  o 

where  second  order  and  higher  terms  have  been  neglected.   Using  these 
relations,  and  the  defining  equation  for  C   ,  we  can  reduce  (9)  to 
the  following: 

Act  =  {-  -^-  sin  a  -  ^-   [C   n  +   CT  ,  6  ]  +  *-  cos(6   -  a  )} 
mv       o   mv    L0    L6   e    v       o    o 
o  o  e       o 

+  {  2—   sin(6   -  a  ) cos  a  -  ■3 —  CT   }  Aa 

v       o    o    mv       o   mv   La 
o  o  o    c 

-  {  *-  sin(6   -  a  )}A9  +  {1  -  -|^2  CT  '    }  A*6  +  ^-  C_  .  sgn(Aa  -(a  -  a  )). 
v       o    o  2mv^  L0        mv   LA  c    o 

o  o  o 

(11) 


Now,  using  this  expression  to  replace  Aa   in  (7) ,  the  equation  for  C,  , 
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and   substituting   the   resulting   expression   into    (10),    yields: 


A6   =4^ 


[■- 


c  T 

Cx*n    ~    C-kax     °"       ~    o Cvr'      [ sin    a 

MO  Mo      e        2v  Ma  mv  o 


-  -^      (CTn  +  CT  .    6    )    -  -S-  cos(9      -a)]} 
mv  L0  L6      e  v  o  o 


-   tCM     +  ~-  C„.    [  &-  sin(9   -  a   )    -  ~-     cos   a^    -  ^-  CT  Aa 


Ma        2v       Ma       v 
o  o 


o        o  mv 


o       m       La 
v  c 

o 


+   (fe   CM.    sin(6      -   a    )}    AG    -  -~ 

zv6     Ma  o  o  2v  M0  Ma 

o  o  o 


<c„:  +  c„,  [i-aS£2cL-]}  A6 


qSc 

2mvz  "Ma  "LA 


2  C„.  C, 


gn(Aa  -  (a  -  a  )) 
c    o 


(12) 


Equations  (11)  and  (12)  then  describe  the  model  we  shall  analyze. 

In  practice,  many  of  the  terms  which  appear  in  (11)  and  (12)  are 
negligible.   Calculations,  using  representative  values,  will  show  that 
for  many,  if  not  all  cases,  little  is  lost  by  assuming: 


£Sc     . 

2mv2   L6 
o 


«   1  , 


2v    Ma 
o 


«   1   , 


and  also  that, 


_  as 


■&—  sin(6   -  a  ) cos  a        CT 

v       o    o    mv       o   mv   La 
o  o  o    c 


«   1   , 


and, 


•^  sin(9   -  a  ) 
v       o    o 
o 


<<   1   . 


Also,  we  have  not  yet  used  our  assumption  that  the  aircraft  was  trimmed 
for  flight  at  a   ,6   .   This  assumption  says  that  when 


Act  =  A6  =  AG  =  0  , 

then, 

Act  =  A6  =  0. 

Thus,  imposition  of  the  trim  conditions  reduces  (11)  and  (12)  respectively 

to: 

0  =   {-  —  sin  a     -  ^-   [CT  A  +  CT  x    6    ]   +  &     cos (6      -  a   )    } 
mv  o       mv  L0  L6      e  v  o  o 

o  o  e  o 


-  ^     CT  .    sgn(a     -  a   )  ,  (13) 

mv        LA      6        c  o 

o 


and, 


0  ■   "  CM0  "  Cm    Se        •  <"> 

e 

If  we  now   incorporate    (13) ,    (14) ,    and  our  assumptions   about   small 

terms   into    (11)    and    (12) ,   we  arrive  at 

Ad  =   AG  +  ^—     CT  .    {sgn(Aa  -    (a     -  a   ))   +  sgn(a     -  a   )}    ,  (15) 

mv       LA  co  co 

o 

and, 

c-2 

AG   =   -  Sf£  r        Aa   -  ^V     {C    '   +  CL,.    }    AG 
I  ■      Ma  2v   I  MG  Ma 

y  o  y 


-   oqSTC)2     CM.    CT.    {sgn(Aa  -    (a     -  a   )   +  sgn(a„   -  a   )  }    .    (16) 
2ml  v         Ma     LA  coco 

y  o 


10 


We  now  define  the  quantities   3,  w,  a,  and  z,     by: 


*f£  CM  =  (62  +  W2)  >  0   , 
I    Ma 


f^~  {C  .  +  CM.}  =  26  >  0  (17) 

2v  I    M6    Ma  v   y 

o  y 


^-     CT   =a>0, 
mv    LA 
o 


and, 


qSc2 


Cx„.  =  5  >  0   . 


2v  I   Ma 
o  y 

Since  it  becomes  important  in  our  later  discussion,  note  here  that 
the  quantity: 


-2 

Sie 


(2e  - » ■  H^  v 


o  y 

Thus  our  model  for  the  flight  dynamics  of  an  aircraft,  trimmed  for 
longitudinal  flight  at  an  angle  of  attack  near  the  stall  break  angle, 
becomes: 

Ad  =  AG  +  a{sgn(Aa-(a     -  a   ))   +  sgn(a  -a   ) }    .      (18) 

CO  CO 

2  2 

A9   =  -    (6     +  a)   )Aa  -  2   3   A6   -   £a(sgn(Aa-(a     _  a   ))   +  Sgn(a  -a   ) } .      (19) 

CO  CO 

(Note  that  (18)  and  (19)  describe  a  system  which  is  second  order  in  Aa 
and   A9  .   Since  the  original  system  (1)  -  (4)  was  equivalent  to  a  fourth 
order  system  in  a,  6,  6,  and  v,   our  model  will  obviously  not  describe  the 
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full  range  of  dynamics.   However,  we  will  show  that  it  does  appear  to 
predict  the  known  limit  cycle  behavior.   In  essence,  what  we  are  doing 
in  equations  (18)  and  (19)  is  basing  our  analysis  on  the  influence  of 
the  stall  break  on  the  short  period  mode  of  the  aircraft,  and  neglecting 
the  long  period  mode.   Furthermore,  the  mathematical  advantages  of  dealing 
with  a  second  order  system  are  enormous.) 

3.   FLIGHT  WITH  TRIM  AT  THE  STALL  BREAK  ANGLE 

We  begin  our  investigation  of  the  model  described  by  (18)  and  (19) 

with  the  case  a  =  a   ,  that  is  the  aircraft  is  trimmed  for  flight  at 
co  ° 

the  stall  break.   Then  since  a  -  a  =0,  the  second   sgn  function  in 

co 

both  (18)  and  (19)  disappears.   If  we  let  x  =  A6   and  x  =  Aa  ,  our 


system  reduces  to: 


jL  =  -  23  x  -  (3  +  co  )  x.  -  Ca  sgn  x. 


x„  =      x  +  a  sgn  x2  (20) 


A.   The  Phase  Plane.   Consider  the  analytic  behavior  of  this  system 
in  the   (x..  ,  x~)  phase  plane.   The  line  x~  =  0  divides  the  plane  into 
two  regions,  where  separately  the  equations  are  linear,  constant 
coefficient  and  non-homogeneous.   For  example,  in  x„  >  0  ,  the  system  is: 

2     2 
x-  =  -  26  x  -  (3  +  to  )x  -  £a   , 

x„  =      x  +  a   .  (21) 
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A  key  observation  now  is  that  (20)  is  symmetric,  i.e. 

{  (x1(t),  x2(t))  |  t1  <  t  <  t2  } 

represents  an  arc  of  a  trajectory  if  and  only  if 

{  (-Xl(t),  -x2(t))  |  t   <  t  <  x2  } 

also  does.   Thus  we  need  direct  our  attention  solely  to  the  region 
x~  >  0,  and  draw  any  needed  information  about  x.  <  0   from  symmetry 
arguments.   We  shall  use  the  usual  convention  in  referring  to  the  line 
x_  =  0  as  the  switching  line. 

We  start  by  considering  the  critical  points  of  the  system,  i.e., 
the  points  where  x  =  x„  =  0  .   These  are  given  by: 

Xl  =  X2  =  °   ' 

(23  -  Q 

xi  =  a'     x2  =  "  ~WT17)a      •  ' 


and, 


(23  -  C) 
Xl  =  "a'     X2  =  (32  +  u,*)   a 


(These  critical  points  actually  occur,  since,  as  we  have  noted, 
(23  -  C)  >  0   ) .   The  presence  of  the   sgn   term  in  (20)  implies  that  the 
critical  point  at  (0,0)  will  be  unstable.   However,  it  is  easily  shown 
that  the  homogeneous  part  of  (21)  has  the  characteristic  polynomial 

X2  +  23X  +  (32  +  w2)   , 
whose  roots  must  have  negative  real  parts  according  to  our  definition 
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of   3  and  co   in  (17).   Hence,  the  two  critical  points  not  at  the  origin 
are  locally  asymptotically  stable.   (That  is,  small  perturbations  off 
this  point  decay  back  into  it.)  We  shall  show  shortly  these  are  not 
globally  asymptotically  stable. 

B.   Conditions  for  a  Limit  Cycle.  We  now  consider  the  conditions 
under  which  (20)  can  have  a  limit  cycle  (i.e.  periodic)  solution. 
Obviously  from  the  stability  discussion  above,  any  limit  cycle  trajectory 
must  pass  through  the  switching  line.   Furthermore,  it  must  by  symmetric 
with  respect  to  the  origin.   Thus,  if  we  consider  the  arc  of  the  limit 
cycle  in  the  region  x~  >  0,  and  assume  this  arc  originates  on  the 
switching  line  at  the  point   (x*,   0),  then,  after  a  half-period,  it 
must  terminate  at  the  point   (-x*,  0).   Thus,  if  we  denote  the  half-period 
by  T  ,  we  have  a  necessary  condition  for  a  limit  cycle  of  period  2T  is: 


x1(0)  =  -  xx(T)   , 


x2(0)  =  x2(T)  =  0  . 


(22) 


(This  situation  is  represented  in  figure  4.) 

It  is  easily  shown  that  in  the  region  x_  >  0  the  solution  to 

(20)  is  given  by: 

— 8t 
x  (t)  s  -c-  e   (3  cos  cot  +  co  sin  cot) 

— 8t 
+  c?  e   (co  cos  cot  -  3  sin  cot)  -  a   ,        (23) 


and, 


x_(t)  =  c.  e    cos  cot  +  c2  e    sin  cot  +  j  2    ■   2\  a  •     (24) 
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x#*  =  switching 
line 

=  critical 
point 


The  Limit  Cycle  Necessary  Condition 
Figure  4 


<8>  =  critical 
point 

mhk   switching 
line 


tXXXXXXXKXXXXXXXKXXXXXXXXXXXXXXXX  fc» 


The  x-  >  0  Region 
Figure  5 
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(These  simply  represent  arcs  of  counterclockwise  spirals  centered  at  the 

(  7  R  —  t  ^ 
point   (-a,  :■  "21  ■  ■ — 2\  a  )  >  which  is  the  critical  point  in  the  upper  half 

plane.   Thus  every  curve  originating  on  the  switching  line  must  follow 
a  trajectory  that  carries  it  above  this  point.   (A  typical  such  curve  is 
labeled  (1)   in  figure  5.)   Furthermore,  since  the  critical  point  (if  it 
exists)  is  asymptotically  stable,  there  must  be  some  curves  entering  the 
region  which  never  exit  it,  but  decay  down  to  the  critical  point.   (One 
of  these  is  shown  as   (2)   in  figure  5.)   Finally,  since  any  trajectory 
intersecting  the  line  x..  =  -a  has  a  horizontal  slope  at  the  point  of 
intersection  (see  eqn  (21)),  then  the  boundary  between  the  decaying  curves 
and  those  that  exit  to  x_  <  0  must  just  be  that  curve  which  terminates 
at   (-a,0)  .   (This  curve  is  labeled   (3)   in  figure  5.) 

But  also  note  that  the  time  interval  represented  by  a  trajectory  in 
the  region  x~  >  0  is  a  function  of  the  counterclockwise  angle  (as 

( 26  —  L ^ 

seen  from   (-a,  -7-77 — i 2\  a  )  )   between  the  initial  point   (x  (0) ,  0), 

{is     +   co  ;  1 

and  the  point  at  which  the  trajectory  first  touches  the  switching  line 

again.   Thus,  any  curve  which  represents  a  half-period  of  a  limit  cycle 

must  spend  less  time  in  x„>  0  than  T   ,   where  T   is  the  time 
r  2  o  o 

represented  by  the  curve  (3)  in  figure  5.   Thus  T  satisfies: 

T  <  T 

where, 

x.  (T  ) '=  -a 

1  o 

x0(0)  =  x_(T  )  =  0  .  (25) 


o 
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Conditions  (22)  and  (25)  can  now  be  seen  to  provide  both  the  necessary 
and  sufficient  conditions  for  existence  of  a  limit  cycle  of  period   2T  . 
The  condition   x  (0)  =  x?(T)  =  0  ,   applied  to  (24),  immediately 


leads  to: 


and, 


(23  -  Q      a  /0„ 

Cl  "  "  WTl^)    a   '  (26) 


RT 
c  (e   -  cos  coT)  =  c     sin  wT  .  (27) 

Notice,  according  to  these,  that   T  ^  —  . 

w 

The  second  condition,  on  x  (t) ,   when  applied  to  (23),  implies: 

x  (0)  =-3c1  +  a)c2-a=-  x  (T) 

— RT 
=  c1e    (6  cos  (joT  +  a)  sin  coT) 

—  RT 
-  c«e    (3  cos  wT  -  co  sin  coT)  +  a 

or, 

—  RT  — RT 

c.  (3  +  3e    cos  coT  +  toe    sin  wT) 

-  c  (oj  +  coe~6T  cos  wT  -  3e    sin  wT)  =  -  2a   .      (28) 
However,  we  have  already  commented  that 

IT 

T  ^  —  (because  (23  -  C)  4   0) .   Then  sin  u)T  ^  0,  and  hence  we  can  multiply 

RT 
both  sides  of  this  last  equation  by  the  term  e   sin  coT,  and  simplify 

{using  (27)  to  replace  the  term  c„  sin  wT)  to  yield: 

/0      3T      .  23TN  „      3T 

c.  (23e        sin  coT  +  to  -  toe        )    =   -   2ae        sin  coT    , 


or, 


(2a  +   &c±) 
sinh   3T  a  sin  uT      , 
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which,  using  (26),  reduces  to: 

Equation  (29)  expresses  the  necessary  condition.   In  order  to  guarantee 
sufficiency,  i.e.,  in  order  that  a  solution  of  (29)  must  generate  a 
limit  cycle,  we  must  add  the  restriction, 


T  <  T    , 
o 


where  T   was  defined  above, 
o 


Applying  the  definition  of  T   to  the  solutions  given  by  (23) 


and  (24)  yields 


(23  -  O 

Cl  "    (6Z+  o)z)  a    » 


cn  (e*  °  -  cos  *jT  )  =  c„  sinuiT    , 

I  O      l  o 


and, 


—  RT 

-cne   °  (B  cos  coT  +   sin  wT  ) 
I  o  o 

—  RT 
+  c„e   °(o)  cos  wT  -  B  sin  wT  )  =  0   . 
2  o  o 

Note  that  the  first  two  equations  are  identical  to  (26)  and  (27)  .   This 
follows  since  the  conditions  on  x„(t)   implied  by  both  (22)  and  (25) 

are  identical.   Again  we  can  assume  sin  wT  #0,  multiply  both  sides  of 

RT 
the  last  equation  by  e  °  sin  coT  ,   simplify  by  raplacing  the  term 

c9  sin  a)T   by  its  value  in  terms  of  the  middle  equation,  and  divide  by 


c, ,   to  yield; 


—  RT 
o)  cos  wT  -  3  sin  wT  =  u  e   o   .  (30) 

o  o 
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In  view  of  our  earlier  comments  then,  we  must  pick  the  smallest 
positive  solution  of  (30). 

Mathematically,  then  we  can  summarize  our  results  in  the  following: 
Theorem:   The  system  (20)  will  have  a  limit  cycle  of  period   2T   if  and 
only  if  T   satisfies: 

sinh  BT  "  "  "(2B  -\~)   "'  Sln  bT>  (29) 


and, 


where: 


T  <  T   , 
o 


T   is  the  smallest  positive  solution  of 
o 

—  RT 

co   cos  coT     -   3   sin  coT     -  to  e        °    .  (30) 

o  o 

C.   Example  -  The  Case  Considered  by  Frederiksen  [1].   To 

determine  validity  of  the  predictions  of  this  model,  consider  the  case 

used  by  Frederiksen  [1],   Figure  6  shows  the  table  of  values  appropriate 

to  this  case,  which  describes  the  F-94,  a  straight  wing,  single  engine 

jet.   The  value  of  the  coefficients  appropriate  to  the  linearized 

differential  equations  were  computed  and  are  shown  in  figure  7.   Note  that 

our  assumptions  on  small  size  of  certain  coefficients  seems  valid.   Using 

the  values  of   3>  co  and  t,      so  computed,  equations  (29)  and  (30)  were 

graphically  and  numerically  solved,  and  are  graphically  portrayed  in 

figures  8  and  9,  respectively.   There  are  obviously  two  solutions  of 

TT  2tT 

(29)  in  the  region  —  <  T  <  —  (which  is  the  region  of  interest  since 

CO  CO 

(2B-C)  >  0.)   These  are 

T  =  1.87      and   T  =  3.56 
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II 

Parameters  for  F-94 


p  =  2.38  x  10  3  c  =  6.40 

m  =  384  I  =  2.65  x  104 


S  =  239.0 


\l     Values  used  by  Frederiksen  [1] . 


y 


CL0  =  1.440  CL6  =  0.43 

e 

CT    =0.30  C*  =  10"5 

La  L6 

c 


CTA  =  0.05 

LA 


SlO  '   °-15  CMa  "  4"27 


Sfa-1"54  CM«  =0-88 


e 


Sie  =  8-16 


9  =  0.099  rad  v  =  169  ft/sec 

o  o 

a  =  0.41   rad  T  =  0.0 

o 

a  =0.41   rad 
c 


Figure  6 
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Based  on  the  parameters  in  figure  6, 

qSc  - 

7:     C*      =   2.4  x  10    ,   -£-  C  .   =   0.081 

0    I        L6  Ma 

2mv  2v 

o  o 


^-  sin  (6   -  a  )  =  -  .058 
v         o    o 
o 


-£-  sin  (0   -  a  ) —  cos  a   -  ^ —  r  =   _  .096 

v        o    o    mv         o   mv   La 
o  o  o    c 


qSc 

- —  CM   =  (3  +  a>  )  =  3.00 

I    Ma 

y 

s-2 

q  °    {C  •  +  C  .}  =26=  0.46,  so  6  =  0.23,  w  =  1.72 

M6    Ma 


2v  I 
o  y 


qs  _3 

C,  .  =  a  =  6.3  x  10 


mv    LA 
o 


s72 


C„.  =  C  =  0.16 


2v  I    m 
o  y 

(2B-0  =  0.30  >  0. 


Coefficients  for  the  Differential  Equations 


Figure  7 
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Figure  8 


Graphical  Solution  of  (30) 
Figure  9 
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But  T  ,   defined  as  the  first  positive  intersection  of  the  curves 
o 

in  figure  9,  is  given  by 

T  =3.10 
o 

Thus,  only  the  lower  value  of  T   satisfies  the  requirement 

77 

—  <  T  <  T   and  there  will  be  exactly  one  limit  cycle,  with  period 

to        o 

2T  =  3.74  sec 

Computing   c.   and   c„   according  to  (26)  and  (27)  ,  and  using  the 

values  to  describe  x_(t)  =  Aa(t)   as  given  by  (24),  leads  to  an 

amplitude  for  this  limit  cycle  of 

Act    =  .018  rad  =  1° 
max 

These  values  are  in  very  good  agreement  with  those  obtained  by 

Frederiksen,  which  were   2T  =  3.77  sec 

Act    =  .020  rad   . 
max 

D.   Additional  Comment.   We  shall  make  one  other  fairly  important 

observation  on  the  mathematical  nature  of  (29)  and  (30) .   Note  that  (29) 

is  equivalent  to 

sinh   3T  qj2  +   SC   -    32  /3n 

sin  wT  '    0>(2B-O 

Now  since   (28-<;)  >  0,  the  region  of  interest  for  this  equation  is 

TT  271 

—  <   t   <  —     .      Letting 

CO  U) 

G(t)    =   Sinh    gt      ,  (32) 
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we  have 


and, 


G(t)    <    0,    in     -  <   T    <   —        , 

CO  CO 


^  /IT    +■,  „  /2lT-. 

G<w  >  -  G<- > "  -  - 

Furthermore 

„w    v         B   cosh   Bt   sin  cot   -  co   sinh   Bt   cos  oat 
«   v.t;   -  2 

sin      cot 

Note  the  presence  of  the  square  in  the  denominator  implies  that  the 
algebraic  sign  of   G'(t)   is  the  same  as  the  sign  of  the  numerator. 
This  numerator, 

B  cosh  Bt  sin  cot  -  co  sinBt  cos  cot 


has  the  value 


and  the  value 


.   i   BfT     n  TT 

co  sinh  —  >  0  at  t  =  —  , 

CO  CO 


.  .  2Btt   n  2tt 

-co  sinh  <  0  at  t  =  — 

CO  CO 

The  derivative  of  this  numerator  is 

2    2 
(B  +  co  )  sinh  Bt  sin  cot. 

,,,..,  .  TT  2lT      ,      .     .,     - 

Since  this  derivative  is  negative  on  —  <  t  <  —  ,  then  it  clearly 

b  CO  CO  J 

TT 

follows  that  the  numerator  of  G'(T)   starts  positive  at  —  ,   remains 

r  co 

positive  up  to  some  point  T   ,  —  <  T   <  —  ,  then  turns  negative, 

r         r         r       m   co    m    co 

staying  negative  for  the  remainder  of  the  interval.   Based  on  earlier 
comments,   G'(t)   behaves  the  same.   Thus  G(t)   is  monotonically 
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7T 

increasing  on  —  <  T  <  T   ,  and  monotonically  decreasing  on 

o  o 

T   <  T  <  —  .   Hence   G(t)   has  a  unique  maximum  in  —  <  t  <  — 
m       co  a)       to 

[Figure  10  shows  the  general  form  of   G(t).] 

But  now  let's  investigate  the  relative  positions  of  T   and  T 

o  m 

At     T    ,    G'(T   )      has   the  numerator 
o  o 

3   cosh   3T      sin  coT     -  co  sinh   3T     cos  coT 
o  o  o  o 

I  Qrn  -|  Qrri 

=  —  e     o    [3   sin  coT     -  co   cos  coT    ]   +  -  e        °    [  3   sin  coT     +  co  cos  coT    ]  . 
I  o  o  I  o  o 


But   according   to    (30) , 


—  RT 
3   sin   coT      =   -   coe        o  +  co   cos   coT 
o  o 

Thus,  the  numerator  of   G'(T  )   is 

o 

1      3Tn    r        -3T    ■,        1      -3T  -  T0        0  •, 

=  —  e     °    [-coe        o]    +  -r  e        o    [-coe  +  2cocos   coT    J 

.  .        —  RT 

=  -r  e        o    [cos   coT      -   cosh  coT    ]    <    0      . 
Z  o  o 

Then 

G'(T  )  <  0    ,     i.e.,   G(T)   is  decreasing  at  T   ,   and  so  by 

our  earlier  comment,   T  <  T  . 

m    o 

But  observe  that  the  right  hand  side  of  (31)  is  independent  of  T  . 
Thus,  since  (29)  and  (31)  are  equivalent,  a  necessary  and  sufficient 
condition  for  (29)  to  have  a  solution,  based  on  our  analysis  of  G(T)  ,  is 


["sinh  Btl 
tt   Lsin  wt  J 


2  2 

u)     +  3^  ~  3 

max  — : —         >     *      > —        .  (33) 

2„       |  sin  cot  00(23-0 

-  <  -t  <  ~ 

0J  CO 
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Figure   10 
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Figure    11 
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Finally,  since  G(-  )  =  -  «  and  T   <  T  ,  then  if  (33)  is  satisfied, 

o)  mo 

then   there  must  be   a     T   <   T     <   T        such   that 

m  o 

2  2 

sinh   BT  _    _  03     +   Be   -   B 

sin  a)T  o)(2B-c) 

Then  this  T   satisfies  (29)  and  (30),  and  so  generates  a  limit  cycle. 
Thus,  (33)  express  a  necessary  and  sufficient  condition  for  existence  of 
a  limit  cycle  in  this  model. 

This  is  expressed  by  the  following: 
Theorem:   The  system  (20) ,  where   (2g-j;)  >  0  ,   will  have  a  limit  cycle 
solution  if  and  only  if 

u2  +  Be  -  B2 
max  — : ■ —    >   -  - 


[  sinh  Bt 


oo(2B-0 
IT        2tt 

—  <  t  <  — 

Observe  that,  with  the  substitution 


X  =  -  >  0  ,   and   n=~>0, 

03  03 


the  necessary  and  sufficient  condition  becomes 
max 

TT<t<2JT 


sinh  Xt      1  +  nX  -  X  ,„■> 

[sin  t   J  "     (2X  -  n)      ' 


(the  left  hand  side  of  which  is  clearly  a  function  of   X   only.)   Note 
we  already  have  commented  that   (2B  -  O  >  0,  hence 

(2X  -  n)  >  0,   or  0  <  n  <  2X   . 
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Now,  observe  that 


d_ 

dn 


1  +  nA  -   A2 

[    (2X  -  n)  J 


a2  +  i 
(2A  -  n) 


2    <  0 


and  therefore  the  right-hand  side  of  (3A)  has  its  maximum  value  at 

C  =  0.   Thus,  clearly,  a  sufficient  condition  for  (34)  to  be  satisfied, 

and  hence  for  a  limit  cycle  to  exist,  is, 


max 

7T<t<2TT 


sinh  At 


sint 


X2-l 

2A 


In  figure  (11)  we  plot  the  function 


max 
TT<t<2ir 


sinh  At 


sint 


1  -  A' 
2A 


for  0  <  A  <  1,   and  observe  that  (35)  cannot  be  satisfied  if 

A  >  .280  . 

(Observe  (35)  cannot  be  satisfied  for  A  >  1.)   Thus,  a  sufficient 
condition  for  existence  of  a  limit  cycle  is: 

A  <  .280 

Now,  if   A<   .280  ,   and   3  =  Xw  , 

we  have: 

2    2    2         2    2  2 

03  <  a  +  3  =  (1  +  X  )  u  <  1.078  a)   , 


fl  2 

and  so  we  can  approximate,   oj  =  v  w  +  B 

o 

But  then,  according  to  the  definitions  (17),   A .=  —  is  proportional 


(35) 


(36) 


to 


CMa  +  CM0 


JMa 


28 


and  we  see  our  sufficient  condition  (35)  translates  into  a  statement 
that 

C    +  C  " 

Md    M6 

is  somehow  small  compared  to    C.. 

r  Ma 

Since  C.,  is  analogous  to  a  spring  constant,  and  C_  *   and  Cw.   are 

Ma         °  r  M6        Ma 

damping  terms,  we  can  view  (35)  and  (36)  as  saying: 

"A  sufficient  condition  for  existence  of  a  limit  cycle  behavior 

in  the  presence  of  a  sharp  negative  lift  slope,  is  that  the  damping 

terms  in  C_.  are  small  relative  to  the  "restoring"  term,   Cw  ." 
M  Ma 

Clearly  this  interpretation  makes  physical  "sense". 

4.   SUMMARY  AND  CONCLUSIONS 

In  this  paper,  we  have  developed  and  analyzed  a  simple  two-dimensional 
model  for  post-stall  time  histories  in  the  presence  of  a  sharp  lift  loss 
at  stall.   This  model  yields  a  simple  necessary  and  sufficient  condition 
for  the  occurrence  of  limit  cycle,  i.e.,  "rocking  chair"  stall.   Further- 
more, the  model  predicts  values  that  are  in  very  close  agreement  with 
those  obtained  from  numerical  solution  of  the  full,  four-dimensional,  non- 
linear case.   Lastly,  the  model  yields  physically  satisfying  information 
about  the  important  parameters  for  occurrence  of  this  limit  cycle. 
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APPENDIX.   Definition  of  Variables  and  Coefficients 

c  =  Mean  Aerodynamic  Chord  (MAC) 

C  =  Drag  Coefficient 

C  =  Lift  Coefficient 

Li 

C„,  =  Pitching  Moment 

M 

g  =  Gravitation  Acceleration 

I  =  Moment  of  Inertia  about  e.g. 

y 

m  =  Mass 

,     2 

q  =  H.   P  v 

5  =  Wing  Reference  Area 
T  =  Thrust 

u  =  Velocity  Component  along  Thrust  Line 

v  =  Velocity 

w  =  Velocity  Compenent  normal  to  Thrust  Line 

a  =  Angle  of  Attack 

6  =  Pitch  Angle 

p  =  Density  of  Air 
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